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Objective Bayesian Inference for Bilateral Data 


Cyr Emile M’lan* and Ming-Hui Chen^ 


Abstract. This paper presents three objective Bayesian methods for analyzing 
bilateral data under Dallal’s model and the saturated model. Three parameters 
are of interest, namely, the risk difference, the risk ratio, and the odds ratio. 
We derive Jeffreys’ prior and Bernardo’s reference prior associated with the three 
parameters that characterize Dallal’s model. We derive the functional forms of the 
posterior distributions of the risk difference and the risk ratio and discuss how to 
sample from their posterior distributions. We demonstrate the use of the proposed 
methodology with two real data examples. We also investigate small, moderate, 
and large sample properties of the proposed methodology and the frequentist 
counterpart via simulations. 

Keywords: Bayes factor, Dallal’s model, Jeffreys’ prior. Odds ratio, Product tri¬ 
nomial distribution, Reference prior, Risk difference, Risk ratio. 


1 Introduction 

Bilateral data arise in medicine when a group of randomly chosen patients with a con¬ 
dition receive a new treatment, for example, surgery, on paired body parts within the 
same individual (eyes, ears, breasts, arms, hands, knees, legs, or feet), while another 
group of patients with this condition receive a control treatment, for example, the cur¬ 
rently accepted medical treatment. The investigator records paired Bernoulli outcomes 
about a particular characteristic, for example, absence of the condition, that are then 
grouped into one of three categories, (i) The two body parts are cured, recorded as (1,1). 
The counts of patients with this characteristic from the control and treatment groups 
are denoted by TO 20 and TO 21 - (h) One of the two body parts is cured while the other 
remains diseased, recorded as (1,0) or (0,1). However, these (1,0) and (0,1) outcomes 
are thrown away and this detailed information is no longer available. Only the counts 
of patients from the control and treatment groups in that category, mio and mu are 
available, (hi) Neither of the two body parts are cured, denoted by (0,0). The counts 
of patients in this category from the control and treatment groups are denoted by moo 
and moi- The data, denoted by ZJ, can be summarized into a 3 x 2 contingency table 
(see Table 1) where the trinomial counts (moi,mu,m 2 i) for the treatment group and 
(moo, mio, m 2 o) for the control group are the cell entries. Such data are very common 
in ophthalmologic, orthopaedic and otolaryngologic studies. Twin studies are also a fa¬ 
miliar source of bilateral data. The goal of such clinical trials is to quantify the benefit 
of treatment over placebo. 
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Table 1: Data structure for bilateral data along with their corresponding trinomial 
probabilities, poo + Pio + P20 = 1 = Poi + Pii + P21, moo + ™io + W 20 = m+o and 
TOoi + mil + m 2 i = m_|_i are fixed by design. 


Numbers of 
cured organs 

Group 

Treatment 

Control 

0 

moi (poi) 

moo (Poo) 

1 

mil (pii) 

mio (pio) 

2 

TO 21 (P2l) 

m20 (P20) 

Total 

m+i 

m+o 


The main parameter in bilateral models is the difference between the proportion 
of body parts (Aq, Ai) with the characteristic of interest in the two groups, A = Ai — 
Aq. However, the dependency between paired observations cannot be ignored. Rosner 
(1982), Morris (1993), and Tang et al. (2006) discuss the consequences of ignoring this 
dependency. 

Five models have been proposed for bilateral data: Rosner’s model, Dallal’s model, 
the equal correlation model, the independence model, and the full or saturated model. 
Of these five models, the most extensively studied is Rosner’s model. Tang et al. (2008) 
present several test statistics for the equality of Aq and Ai under Rosner’s model. Qiu 
et al. (2009) consider the problem of sample size calculations under Rosner’s model. Tang 
et al. (2011) discuss various techniques to construct asymptotic conhdence intervals 
for A under Rosner’s model and evaluate the performance of these via empirical studies. 
Tang et al. (2010) is perhaps the only paper that discusses all five models simultaneously 
to determine which one provides a better ht to the data. Pei et al. (2010) present 
asymptotic confidence intervals under the equal correlation model and evaluate their 
performance. Pei et al. (2008) present a test for the equality of Aq and Ai when unilateral 
and bilateral data are combined under the equal correlation model. 

Among these five models, the independence model is rarely used in practice in the 
context of bilateral data. The objective priors for the saturated model have already been 
investigated in the literature. The challenge with Rosner’s model is that it is difficult to 
justify this model from a biological point of view. Only the equal correlation model and 
the equal conditional probability model (Dallal’s model) have sound statistical foun¬ 
dations and biological interpretation. Indeed, one can always characterize a bivariate 
discrete distribution for two binary random variables by making assumptions about (i) 
their marginal distributions and the correlation they share or (ii) their marginal distri¬ 
butions and the two conditional distributions they share. For Rosner’s model and the 
equal correlation model, Jeffreys’ prior and Bernardo’s reference prior are too complex 
to be of much practical use. Dallal’s model is the only model for which we could derive 
useful closed-form expressions for various objective priors. For these reasons, we solely 
focus on Dallal’s model in this paper. 
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Bayesian methods have gained incredible popularity in recent years both in the the¬ 
ory and practice of statistics. Under non-informative priors, Bayesian inferences yield 
results similar to that obtained under the frequentist paradigm. As we will show in an 
example presented in Section 5.2, Bayesian inference for Dallal’s model yields similar 
findings to the corresponding frequentist analysis. Bayesian methods, however, do not 
rely on the normal approximation to carry out statistical inferences, which is an advan¬ 
tage over the frequentist methods. It is not unusual to encounter 3x2 bilateral data 
where one or more cells have sparse data, thereby preventing the use of the usual nor¬ 
mal approximation that underlies frequentist inferences developed in Appendix 1 of the 
Supplementary Web Materials. Also, in some sparse 3x2 bilateral data, frequentist es¬ 
timates lie on the boundary of the parameter space, which are not permitted by design. 
In such situations, one cannot compute the confidence intervals for some parameters 
or carry out tests of hypotheses about some of the parameters. We provide one such 
example in Section 5.1. Bayesian methods provide a simpler way to analyze such sparse 
bilateral data. Another benefit of Bayesian inference for 3 x 2 bilateral data is its abil¬ 
ity to handle the nuisance parameter in Dallal’s model which complicates frequentist 
analyses. 

There are no existing Bayesian methods for bilateral data in the literature yet. We 
present objective Bayesian inferences for three parameters of interest: the risk difference, 
the risk ratio, and the odds ratio. In Section 2, we present Dallal’s reduced model along 
with Dallal’s full model. Section 3 is dedicated to the derivation of the objective Bayesian 
modeling of bilateral data. We focus primarily on deriving Jeffreys’ prior and Bernardo’s 
reference prior. We then discuss the Bayes factor in the context of hypothesis testing 
as well as a simulation scheme for the joint posterior distribution. Section 4 presents 
results of an empirical comparison between Bayesian methods and frequentist methods. 
Section 5 presents two illustrative case studies. In Section 6, we present two families of 
Bayesian prior distributions, including Jeffreys’ or Bernardo’s reference priors as special 
cases. Section 7 concludes this paper. 


2 Dallal’s Dependence Model 

Dallal’s model was presented for the first time as model 2 in Dallal (1988). Let “i = 1” 
and “i = 0” denote the treatment group and the control group, respectively. Denote 
by Zijk a binary variable such that Zijk = 1 if the kth site (fcth body part) of the 
jth subject in the zth treatment group is free of disease at the end of the study and 
0 otherwise for f = 0, 1, j = l,...,m+i, and k = 1,2. Dallal’s (reduced) model is 
characterized by the following assumptions: 

Assumption 1: Pi^Zijk = l) = Ai for z = 0, 1 with 0 < Ai < 1. 

Assumption 2: P(^Zijk = 1 | = l) = 1 — 7 with 0 < 7 < 1. 

Assumption 2 states that the conditional probability of an occurrence of a particular 
characteristic at one site given an occurrence of that characteristic at the other site to 


142 


Objective Bayesian Inference for Bilateral Data 


be the same in the two treatment groups. This statement is relaxed and replaced by 
P(^Zijk = 1 I = l) = 1 — 7 i with 0 < 7 i < 1 , i = 0,1 in the full or saturated 

model. That is, two conditional probability statements are made, one for the treatment 
group and the other for the control group. We also refer to this saturated model as 
Dallal’s saturated model. However, the full model has one more parameter than the 
reduced model. 

Let ruhi be the number of subjects in the ith group with exactly h site(s) cured and 
Phi be the success probability associated with ruhi for h = 0,1, 2 and i = 0,1. The two 
group total sample sizes are denoted by m+i for the treatment group and m+o for the 
control group and these are assumed fixed by design. Hence, (moi, mu^ m, 2 i) follows the 
trinomial distribution with total number of trials m+i and probability parameter vector 
{PQi,Pu,P 2 i) for i = 0,1 such as 


Poj = 1 - (1 + 7)^i: Pii = 2jXi, and p 2 i = (1 - 'y)Xi- 

Dallal’s model also implies that the correlation coefficients between the Zijk variables 
take the form: 


Pi ~ =1 “ z = 0,1. (1) 

1 Ai 

In fact, the conditional probability assumption in Dallal’s model can be replaced by 
the statement in (1) about the correlation coefficient. The correlation coefficient, pi, 
takes both positive and negative values over the entire range (—1,1). The excess risk is 
defined as 


Si = P{Zijk = 1 I Zij(^s_h) = 1) — P{Zijk = 1) = 1 — 7 — Ai 


for z = 0 , 1 . 

The main parameter of interest in this investigation is the risk difference, A = 

Ai — Ao, and, therefore, 7 can be viewed as a nuisance parameter. The risk ratio, i? = —!- 

Ao 

and the odds ratio, ip = interest. Another parameter of interest 

(1 — Ai)Ao 

is the difference of excess risks in both the treatment and the control groups, Jq — (^i, 
which is equal to A under Dallal’s reduced model. To date, the risk ratio, the odds 
ratio and the difference of excess risks have never been discussed in the bilateral data 
literature. The first two parameters add another dimension to the utility of bilateral 
data so that they can be collected under either a prospective study paradigm or a 
retrospective observational study paradigm, allowing for more applications than those 
under the current clinical setting. Although the addition of these new parameters poses 
more challenges for carrying out frequentist inference, no additional work is required in 
the Bayesian framework. 
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What makes the inferential process challenging in Dallal’s model is that one deals 
with a constrained parameter space. Indeed, the parameter space is 

fl = |( 7 ,Ao, Ai) : 0 < 7 < 1 if 0 < max(Ao, Ai) < 


0 < 7 < - 7T—TT - 1 if max(Ao, Ai) > ^ 

max(Ao,Ai) 2 


or equivalently 


n = |(7, Ao, Ai) : 0 < 7 < 1 and 0 < Aq, Ai < 


( 2 ) 


We adopt the second representation in the sequel. 
The likelihood function can be expressed as 
1 




i=0 


m+i 


[l-(l + 7)A.]™“‘(27A,r-[(l-7)Ay 


—-y)™20+™2i;^mio+m20^mil+m2i + 7) Aq] 

x[ 1-(1 + 7)Ai]’"“\ (7,Ao,Ai)Ga 

In Appendix I of the Supplementary Web Appendix, we develop for the first time large- 
sample frequentist inferences for the risk difference, the risk ratio and the odds ratio 
under Dallal’s model. 


3 Bayesian Analysis 

A key component in Bayesian analysis is the choice of the prior distribution. Tradition¬ 
ally, Bayesians have turned to conjugate priors. However, the concept of conjugate priors 
is only universal in standard univariate problems. In addition, when little or no prior 
information is available, conjugate priors become subjective. Thus, non-informative or 
objective priors are more widely accepted. The uniform distribution over the parame¬ 
ter space is an obvious non-informative prior. Jeffreys’ prior and Bernardo’s reference 
prior are alternative choices that are invariant under any parameterization or a larger 
class of parameterizations. Another alternative choice is to use a joint prior that is the 
compromise between an informative prior and a non-informative prior (Sun and Berger, 
1998 ). 

We discuss four types of priors: the uniform prior, Jeffreys’ prior, Bernardo’s prior, 
and Sun and Berger’s reference prior in light of partial information. We derive the 
posterior distributions of A and R as well as P(A > Aq | D) and P[R > Rq\D) in 
Appendix 4. Although we do not provide the posterior distribution of ip (more complex), 
we discuss in Section 3.4 how to sample from the posterior distributions of A,R and 
Ip. We use the parameter values generated to compute posterior probabilities such as 
P(A > Aq \D),P[R > Ro\D), P{ip > iPq\D'), and Bayesian credible intervals. The 
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Bayes factors are introduced in Section 3.5 to compare the models under the hypotheses 
(i) Ho'. Ao = Ai and Hi: Ai ^ Aq and (ii) ifg- 7 i = 7 o versus 71 jo- 


3.1 The Uniform Prior Distribution 

The prior 7r[/(^, Aq, Ai) = 2 for ( 7 , Aq, Ai) € fl refers to the uniform distribution under 
Dallal’s model. This prior is proper and expresses a complete indifference of one vector 
of parameter values over another. The resulting posterior distribution is 


7 r( 7 ( 7 ,Ao,Ai|iA) oc 2 ™i+ 7 ™i +(1 - 7 )™^+[l - (1 + 7 )Ao] 


x[l- (1 + 7 )Ai]’ 


(7, Ao, Ai) G fl, 


(3) 


where mi+ = mio + tou and TO 2 + = TO 20 + m, 2 i. As a result, the marginal posterior 
distribution of the nuisance parameter 7 is 


2mi+ + l 


B(mi++ l,m 2 ++ 1 ) ( 1 + 7 )'"i++™ 2 ++ 2 ’ 0 < 7 < 1 , 

and the conditional posterior distribution of A^, i = 0,1 given 7 is 


nu(^^|7,^) = (1+7)™“+™-+!4"^ (1 + 7)A*] 


B(mu + m 2 i + 1, moi + 1) ’ 
1-7 


0 < Ai < 


(4) 


(5) 


where B(.,.) refers to the Beta function. In other words, 


1 + 7 


Be(TO2+ +1, mi+ +1) 


and (1 + 7 )Ai |7 ~ Be(mii + m 2 i + l,m-oi + 1), i = 0,1, where the notation Be (a,/?) 
represents the standard Beta distribution with shape parameters a and /3. The uniform 
prior can be viewed as a process of adding 1/2 to the summary statistics in the bottom 
four cells of the 3x2 table and 1 to the top two cells. The uniform prior is appealing 
in situations where the physical system imposes a natural parameterization with a 
nice physical interpretation. In general, the uniform distribution lacks the property of 
parameterization invariance. 


3.2 Jeffreys’ Prior 

Jeffreys’ prior has the property of being invariant under a one-to-one reparameterization 
(Jeffreys, 1946). Regardless of the parameterization used, Jeffreys’ prior distribution is 
proportional to the square root of the absolute value of the determinant of the Fisher’s 
information matrix. 


Define U = (1 + 7 )Ao or equivalently Aq = 


U 


1 


alently Ai = 


V 


■7 


■, and let 1 + = (1 + 7 )Ai or equiv- 


1 + 7 


. Under this new parametrization, the parameter space reduces to 


the interval (0,1) for each of the three parameters 7 , U and V. Thus, we are no longer 
dealing with a constrained parameter space. Moreover, the triplet ( 7 , U, V) forms a set 
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of orthogonal parameters in the sense of Cox and Reid (1987), that is, the off-diagonal 
elements of the expected Fisher information matrix are all zero. Propositions 3.1 and 
3.2 give Jeffreys’ priors for ( 7 , U, V) and ( 7 , Aq, Ai), which are derived in Appendix 1. 

Proposition 3.1. Jeffreys’ prior under the parameterization ( 7 , U, V) is 
7 rj( 7 , u, v) oc — -- (u + ( 6 ) 

(1 +7) 


for 0 < y,u,v < 1 and it is proper, where r 
the two treatment groups. 


m+i 

m+o 


is the ratio of the sample sizes in 


Under Jeffreys’ prior, the nuisance parameter, 7 , is independent of both U and V 
and its marginal prior distribution is given by 


71-J ( 7 ) 


TT ( 1 + 7 ) 


0 < 7 < 1. 


(7) 


However, Jeffreys’ prior depends indirectly on the sample sizes in both groups through 
their ratio, r. 

Proposition 3.2. In the original space, Jeffreys’ prior reduces to 


7 rj( 7 iAo,Ai) oc 


(1 -f 7 )(Ao -f rAi) 


y 7(1 - 7 )AoAi [1 - (1 -f 7 )Ao] [1 - (1 -f 7 )Ai] 

( 7 i-^Oi Ai) £ H. 


( 8 ) 


The posterior distribution resulting from the use of Jeffreys’ prior is 


T^jilAoAilD) oc 


X (1 + 7)Ai] 


30 - 1 - 1 / 2-1 

(7, Ao, Ai) G fl. 


B(toi+ -I- 1/2,7712+ + 1/2) 

X (Ao + rAi)i/2 ^^,o+m,o+i/2-i [1 _ (1 + ^)Ao]™™+'/'-' 

moi-l-l/2 —1 


(9) 


3.3 Reference Priors 

Despite its success in the one-parameter context, Jeffreys’ non-informative prior method¬ 
ology often runs into serious difficulties in multiparameter problems (Datta and Ghosh, 
1996). The prior distribution may be difficult to derive and too complex to be easily 
interpretable. This is often the case when several nuisance parameters are present. For 
our problem, Jeffreys’ prior is difficult to interpret given that it depends on the ratio of 
sample sizes m+i and m+Q. In this context, the reference prior may be more preferred 
(Bernardo, 1981; Berger and Bernardo, 1989), loosely defined as vague priors with the 
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least amount of information. Here, we divide the vector of parameters into two sets: 
parameters of interest and nuisance parameters. Then, we consider the parameters se¬ 
quentially in the process of deriving a reference prior. Berger and Bernardo (1992a,b,c) 
took this idea to another level by suggesting to split the parameter vector into multiple 
groups according to their orders of inferential importance. The two authors provided 
a general algorithm for the construction of reference priors. Hence, reference priors are 
not uniquely defined. We derive the reference prior for the parameters 7 , C/, and H, 
leading to an induced reference prior for 7 , Aq, and Ai. 

We start with the two group orderings: (i) {U,V} and then { 7 } and (ii) {V, U} and 
then { 7 }. 

Proposition 3.3. Bernardo’s reference prior corresponding to the groups ordering 
{U,V} and then { 7 } or {V,U} and then { 7 } is 


2l/2y/2-l(i _ ^)l/2-l yl/2-l(l _ y)l/2-l „l/2-l(l _ ^)l/2-l 

B(1/2,1/2)(1 + 7 ) B(1/2,1/2) B(1/2,1/2) 


0 < 'y,u,v < 1. That is, ^ ^ ^ ^ Be(l/2,l/2), U ^ Be(l/2,l/2) and V ~ 
Be(l/2,l/2). 

Proposition 3.4. In the original parameterization, ( 7 ,Ao,Ai), Bernardo’s reference 
prior is equivalent to 


7r^(7 Ao Ai) = — 71 / 2 - 1(1 - 7)1/2-i^^o^_ 

^7 U 71 B(l/2,l/2) 


1 / 2-1 




B(l/2,l/2) 


( 7 i '^Oi -^i) S H. ( 11 ) 


The proofs of Propositions 3.3 and 3.4 are given in Appendix 2. 

Unlike Jeffreys’ prior, 7 , U, and V are independent under Bernardo’s reference prior. 
The posterior distribution resulting from the use of Bernardo’s reference prior is 


7rfl(7> Ao,Ai|£)) 

„ _ ^ +~*~ / _ mi+-|-l/2-l/i _ \m2+-|-l/2-l\mio-l-"A20-l-l/2-l 

B(mi+ + l/ 2 ,m 2 + + l/ 2 )^ 

X [1 - (1 + ^);^^]™oo + l/2-l^™77+m2l+l/2-l [1 _ (1 + + 

( 7 ,Ao,Ai) € H or equivalently 


7Tii{'y,u,v\D) 

2mi+ + l/2 ..ymi+-|-l/2-l(2^ _ .^^m2+-|-l/2-l 

B(toi+ 1 / 2 , m 2 + -f 1 / 2 ) (1 ^)mi + -Hm 2 + + l 

.ymio+m2o-l-l/2-l(l _ .j^^moo + 1/2-1 .yinii-|-m2l-l-l/2-l (1 _ ^'^moi + 1/2-1 

B(mio-f 17120 + 1/2, Togo- f 1/2) B(toii-I- 77121 -f 1/2, moi-I-1/2) ’ 









Cyr Emile M’lan and Ming-Hui Chen 


147 


0 < 7 ,'u,'u < 1. The reference prior can be viewed as adding 1/4 to each of the bottom 
four cells of the 3x2 table and 1/2 to the top two cells. 

Ghosh and Mukerjee (1992) advise reversing the role of parameters of interest and 
nuisance parameters to obtain a reverse reference prior. That is, reconsider the group 
ordering of { 7 } and then {U,V} or { 7 } and then {V^U}. The reference prior remains 
unchanged as shown in Appendix 2. We also discuss the idea of reference priors under 
the partial information introduced by Sun and Berger (1998) and known as conditional 
reference priors in Appendix 2. 


3.4 Sampling from the Posterior Distribution 

With analytical solutions difficult to derive or compute, we turn to Monte-Carlo simu¬ 
lation methods. We first discuss how to simulate 7 from the distribution 


/(t) = 




1 ' 


M-l (^ _ 


( 1 - 7 ) 


v-1 


B(Ai,i/) (1-^7)''+^' 


0 < 7 < 1 , 


which includes the marginal posterior distributions, t:u{i\D) = 'Kj(yj\D), as special 
cases with /r = mi+ -I- ^ and v = m 2 + + \ - Then, we discuss how to simulate (Ai, Aq) 
jointly from /(Ai,Ao | D). 

We propose two direct and efficient approaches to generate 7 from f{y). 


(i) Let 7 = 


l + e^ 


. We show in Appendix 3 that </ = T — log(2) has the same 


distribution as logit (p) = log 


1 -p 


, where p ^ Be(p, v). Thus, simulate pi 


6 

Be(/i, i/), compute (j)i = logit(pj) — log(2), and set 7^ =- z = 1,..., M. 

1 + 


1 — TT 


(ii) Let 7 = --. We show in Appendix 3 that tt ^ Be(i/,//). Thus, simulate tt^ 

1 + TT 


+ 

■ 

Be(z/, p) and set 7 ^ =- i = 1,..., M. 

1 + TTj 


We now focus on the joint marginal posterior distribution of (U, V) obtained under 
Jeffreys’ prior 


f(u, v\D) (X {u + 


ymio-|-m20-l/2^2 — y)™00-l/2 

B(toio -I- m 20 + 1/2, moo + 1/2) 


B(mii -I- m 2 i -I- 1/2, moi -I-1/2) ’ 


which is independent of 7 . To simulate M observations {ui,Vi), i = 1, • • • ,M, we 
proceed as follows: 
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(a) Simulate independent observations {ui,Vi), i = 1, • • • ,M, with Ui ^ Be(mio + 
77120 + 1 / 2 , Woo + 1 / 2 ) and ~ Be(mii + m 2 i + 1 / 2 ,mpi + 1 / 2 ). 

(b) Compute the weights Wi = {ya + j = 1, • • • , M. 

(c) Use the acceptance/rejection sampling method: Simulate ^ [/(0,1) and accept 
the pair (u^, Vi) only if < Wi/{\ + r). 

(d) Or use the importance sampling method, where all the pairs {ui,Vi) are accepted 
and use the weights Wi to correct for the bias in the computation of posterior 
mean and quantiles. 


Under the reference prior, we simulate independent observations {ui, Vi), i = 1, - ■ ■ , M, 
with Ui ~ Be(mio + W20 + 1/2,moo + 1/2) and Vi ~ Be(mii +m 2 i + 1/2,moi + 1 / 2 ). 
Having simulated a triplet {^i,Ui, U), we compute Aq = Ui/{l + '-fi) and = U/(l + 7 i) 

as well as the risk difference, Ai = ^the risk ratio, Ri = ^ (does not depend 

(l + 7 i) Ui 


on 7 ), and the odds ratio, ipi = 


U[l+7^-tl^] 
U. [1 + 7i - Vi] 


. These simulated values are in turn 


used to compute posterior probabilities and Bayesian credible intervals such as equal¬ 
tailed intervals and highest posterior density (HPD) intervals. Bayesian HPD intervals 
are the shortest intervals containing the parameter of interest with the desired posterior 
coverage probability. They are more desirable than the commonly used equal-tailed 
intervals when the posterior distribution is highly skewed, but are more difficult to 
compute. Chen and Shao (1999) develop the Monte Carlo method to compute the 
HPD intervals. In their paper, they also discussed how to compute Monte-Carlo-based 
Bayesian credible intervals under importance sampling. 


3.5 Marginal Predictive Distribution and Bayes Factor 

In this section, we derive the marginal predictive distribution and the Bayes factor under 
Jeffreys’ prior and Bernardo’s prior. For simplicity and clarity of the presentation, we 
work with the parameterization ( 7 , [/, U). In addition, we derive a single Bayes factor 
formula for these two priors. To accomplish this single formulation, we use a family 
of prior distributions that encapsulates both Jeffreys’ prior and the reference prior as 
special cases. 

3.5.1 Hq-. Xq — Ai versus Hy. Ai ^ Aq 

Consider the hypotheses iJp: Aq = Ai = A and Hi: Ai ^ Aq or equivalently Hq: U = 
V = 9 against Hi: U . Under Hi, we consider the family of prior distributions 

TTffi (7, u, v) 

1 25 '-)/5“l(l — z;3“^(l — 

^KB{^) + Bill) ’ 
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where K is the normalizing constant and 0 < 7, m, r; < 1. Note that when d = 0 , K = 1 . 
Two choices of d are of interest: d = 0 corresponding to the reference prior and d = lj 2 
corresponding to Jeffreys’ prior. The marginal predictive distribution under Hi is 


PHi (rnio, 77120, m- 11 , 77121 ) 


= 77720, TOii, 77121 ) / / (u + rv)'^ 


/o ^0 


jrnio+m2o + ^-l _ y^moo + |-l 

B(?77io + 77120 + f ,77700 + \) 


X —- YT~dvdu, 

B(777ii + 77721 + 2,77701 + 2 ) 


where 

p(777 l0, 77720,77711, 77721 ) = 


777+0 W 777+1 \ B(777i+ + ^,7772+ + \) 

77700,77710,77720/ Woi, 77711,77721 / 


B(ii) 

B(7nio + 77720 + 5,77700 + 5) B(777 ii + 77721 + 5, TTIqi + 5) 


B(i^) 

Under Hq, the likelihood reduces to 

^7,0)=In 


B(i^) 


.2=0 


m+i 

'^2i/ 


2^1+ 


mi+ 


(1-7) 


1712 + 


0mi++m2+ I" 2 ^ _ 01^7170+ ^ 


(1 _(_ j'jmi++m2+ 

where 0 = (1 + 7)A. Jeffreys’ and Bernardo’s priors belong to the family of priors 
2 ^ 75-1 (1 - 7)1-1 (1 - 0)7-1 




B(i,i) (I+ 7 ) B(a,i) 


0 < 7,0 < 1. 


Thus, the resulting marginal predictive distribution is 
PHoimiO, 77720, 77711,77721) 

777+0 \ / 777+1 \ B(777i+ + 7772+ + 

.77100,77110,77720/ V”l-01, TUn , 7772l/ B(i,i) 

B(7772 + + 7771+ + O, 777o+ + |) 

^ BKi) ■ 

The ratio of these two marginal predictive distributions, BFqi, under the condition 
P{Hi) = P{Ho) = ^ (the Bayes factor for testing Hq vs Hi) satisfies 

1 _ B(777io + 77720 + ^, 77700 + ^)B(771ii + 77721 + |, 777oi + ^)B(a, i) 


BF\ B(7772+ + 7771+ + 0,7770+ + i)B( 5, 5)B(i, i) 

f"! 7 I J.y'jd ^17110+77120 + ^-1 _ y^ljlTlOO + ^-l ^17111+77121 + ^-1 _ y^imoi + ^-l 

K B (77710 + 77720 + 5,77700 + 5) B (tTIh + 77721 + ^, 777oi + 5) 


dvdu. 


0 ^0 


Under Jeffreys’ prior, the constant K and the integral term in the Bayes factor are 
computed using computer simulation. Under the reference prior, the integral term dis¬ 
appears and the Bayes factor is computed exactly using only the Beta functions. 
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3.5.2 H*: -fi — 7 o versus H*: yi ^ 70 

One of the statements made in Dallal’s model is that the parameter 7 is constant. As 
discussed earlier, this assumption can be relaxed to = 11 = l) = 1 — 7 ^, 

giving rise to the full or saturated model. Therefore, it is important to test the hypothesis 
Hq : 7 i = 7 o = 7 (Dallal’s reduced model) versus the alternative hypothesis : 73 ^ 79 
(Dallal’s full model). Under H*, U and V are redefined as follows: [/ = (!+ 70 )Aq and 
U = (1 + 7 i)Ai. Under Hq, the prior is TTHi{'y,u,v) defined in Section 3.5.1. Under HJ, 
the family of priors under consideration is 




B(^4)(1+7o) 

ydl-l 

B(ai,A) 


2272 ^(1_7^)2 1 l)l_y)2 1 

B(ii)(l+7i) B(ao,i) 

0 < 7o,7i!'«:'*^ < 1 - 


The resulting marginal predictive distribution is 
Pff* (mio,m2o,mii, 77121 ) = 


B(771io + 7 ,m 2 o + 5 ) 
B(mii + A,to 2 i + I) B(mio + W 20 + ap,mpp + \) 


771+p 

Wip, 77120 


777+1 

77701,77711,77721 


B(i^) 

B(771ii + 77721 + 01,77101 + \) 

B(ai,i) 


B(ao,i) 


We are only concerned with two sets of parameters choices: d = 1/2, oo = ai = 1 used 
in Jeffreys’ prior and d = 0, ao = oi = 1/2 used in the reference prior. 

The Bayes factor for testing vs iJ/ under the condition P{Hi) = P{Ho) = ^ is 


Pffl (777 1 0 , 77720,77711,77721) 

^ 

PH* (mio, m2o, mil, m2ij 
_ B(t77i+ + A, 7772+ + ^)B(ao, A) B(ai, A) 

B(777io + i, 77720 + i)B(777ll + i, 77721 + \) 

^ B(777io + 77720 + ^,77700 + i)B(777ll + 77721 + ^, 777oi + ^) 

B(777io + 77720 + 00,77700 + 5)B(777ii + 77721 + 01,77701 + \) 

I' I 10“1“"2 ^ 


0 -/o 


K 


B(777io + 77720 + 5 , 77700 + 5 ) 

( 1 - 77 ) 


"101 + 4 —1 


B(777ii + 77721 + ^, 777oi + A) 


w} 


dvdu. 


Under Jeffreys’ prior, BF^ is computed using computer simulation while under the 
reference prior it is computed exactly. 















Cyr Emile M’lan and Ming-Hui Chen 


151 


4 Comparisons of Bayesian and Frequentist Intervals: 

An Empirical Study 

In this section, we investigate small, moderate and large-sample performances of fre¬ 
quentist confidence intervals (FCIs) and Bayesian credible intervals (BCIs) under three 
criteria. For a set values for the model parameters, 10,000 3 x 2 bilateral data tables 
are generated from the product of trinomial distributions under a balanced design. The 
essence of these criteria rely on the principle that good FCIs (Wald FCIs described in 
Appendix 1 of the Supplementary Web Materials) or good HPD BCIs should have their 
true coverage close to or preferably larger than the nominal value. Indeed, FCIs and 
BCIs with deflated true coverage are recommended against. Lengths of the intervals 
must also be considered. We use the three following criteria: 

(i) the expected true coverage probability (ETCP) of the interval (A/,, A[/) for A, 
P{Al < A < A[/), in repeated sampling; 

(ii) the expected width (EWCI) of the interval of (Ai,A[/) for A, Au — A^, in 
repeated sampling; and 

(iii) the expected mean square error (MSE) that reflects a compromise between bias 
and precision, in repeated sampling. 

Although these criteria are defined in terms of A, we also examine the behaviors of 
these criteria for ( 7 ,Ao,Ai) as well. For the case of the MSE, we go further by adding 
the three MSEs corresponding to ( 7 , Aq, Ai) to obtain a global measure of MSE. Four 
equal sample size scenarios are chosen: m+o = m+i = to = 10,25,50,100 to re¬ 
flect small, moderate and large sample size situations. For symmetry, we only con¬ 
sider cases where A is non-negative. More specifically, we examine cases with A = 
A/j = h/10, h G {0,1, 2, 3,4, 5, 6 , 7, 8 , 9}. Each choice of A implies the constraint: 
0 < 7 < 7 h,max = min(l, l/A/^-l). We then choose a grid of 7 points: jhj = jjh,max/W 
with j = 1, 2,..., 9. These values are meant to capture a wide range of behaviours of the 
conditional probability of an occurrence of a particular characteristic at one site given an 
occurrence of that characteristic at the other site in the range (1 — 7 ^,max, !)• Then, set 

A?ij,max = - Ah- For the pair (Aq, Ai), we select the values: Xohjk = kXhj,max/W 

i + lhj 

with k = 1,2,... ,9 and Xihjk = Xohjk + A^. Hence, for each sample size and Ah, we 
compute the three criterion functions for 81 combinations of the triplets ( 7 , Aq, Ai). In 
the Bayesian framework, we examine HPD intervals under the uniform, Jeffreys’, and 
reference priors. 

The probability of obtaining degenerate results in the frequentist framework (non- 
estimable model parameters or model parameters on the boundary of the parameter 
space or parameters with a zero variance) is relatively high in smaller samples, and 
even higher when combined with large and small 7 values. These cases are eliminated 
from the frequentist calculations. The five tables in Appendix 2 of the Supplementary 
Web Materials give a summary of results of the simulation study. In each cell, we 
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compute two numbers: the proportion of empirical coverages within 0.01 of the nominal 
coverage and the proportion of empirical coverages above - 0.02 of the nominal coverage. 
The latter summary carries out more value than the former one. The findings for the 
ETCP criterion are summarized as follows. 

(a) For A = 0, Bayesian HPD intervals perform similarly for moderate and large 
sample sizes (m = 50,100) and the proportions of coverages above -0.02 of the 
nominal coverage are close to 100% for Aq, Ai and A. This finding actually holds 
true for A < 0.4. 

(b) Wald FCIs perform poorly for estimating Aq and Ai when m = 10,25 and A < 0.3. 

(c) Overall, in terms of coverage probability, the uniform distribution seems to per¬ 
form better than the other methods when A < 0.5. When A > 0.6, Jeffreys’ prior 
and the reference prior perform better when estimating Aq, Ai and A. 

(d) When it comes to estimating 7 , the uniform prior outperforms the other methods 
regardless of the nominal value of A. 

(e) In general, Jeffreys’ prior and the reference prior tend to give similar results re¬ 
gardless of the nominal value of A. 

A summary of the findings for the MSE criterion is given as follows. 

(a) Jeffreys’ prior and the reference prior perform similarly regardless of the nominal 
value of A. The MSEs for A are scattered around the line y = x, with increasing 
deviations as n gets smaller. As A increases, the deviations from the 45 degree line 
also increase with a slim advantage of Jeffreys’ prior when A > 0.5. See Figures 
1 and 2. The same pattern is observed with the global measure of MSE obtained 
by adding the MSEs for 7 , Aq and Ai. 

(b) Jeffreys’ prior and the uniform prior have equivalent properties for moderate and 
large sample sizes when A < 0.6. However, when A > 0.7, the superiority of 
Jeffreys’ prior is highly evident. 

(c) Wald FCIs perform worst when A < 0.6, with differences worsening as n gets 
smaller. This behavior is reversed when A > 0.7. See Figures I and 2. 

(d) In terms of the MSEs corresponding to 7 , all the Bayesian approaches perform 
better than Wald’s approach when A < 0.8. 

The findings for the EWCI criterion are the following. 

(a) Jeffreys’ prior and the reference prior give similar results when estimating A. See 
Figure 4. The reference prior tends to be better than the uniform prior when 
A > 0.6 and inferior when A < 0.5. The Wald FCI is the worst when A < 0.6 
and the best when A > 0.7. See Figure 3. 

(b) Wald FCIs perform worst when estimating 7 . The uniform prior is no better than 
Jeffreys’ prior with a slight edge to Jeffreys’ prior when the sample size is small. 
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A=0.0 and m=100 A=0.0 and m=50 





Figure 1: Graphs of the 81 empirical MSEs for A generated from the case A = 0. 


5 Case Studies 

5.1 Bayesian Analysis of Bilateral Data with Sparse Data 

Mandel et al. (1982) considered a double-blind randomized clinical trial which compared 
the antibiotics Cefaclor and Amoxicillin for the treatment of otitis media with effusion 
(OME). Among a total of 214 participants in the trial, 11 children were at least six 
years old and underwent bilateral tympanocentesis prior to randomization into one of 
two groups (Cefaclor or Amoxicillin). Children in each treatment group received a 14- 
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A=0.8 and m=50 A=0.8 and m=10 








ms8 for A under Uniform prior rrkse for A imder Uniform prior 

Figure 2: Graphs of the 81 of empirical MSEs for A generated from case A = 0.8. 


day course of treatment with one of the antibiotics and dichotomous ear outcomes were 
determined (i.e., cured or not-cured) and recorded. Table 2 provides a summary of the 
data collected. The primary goal of this investigation is to test if the cure rates were 
identical between Cefaclor and Amoxicillin. A further goal is to estimate the size of the 
difference of the percentage change, i? — 1, in the performance of the two medications. 
This example was discussed in Tang et al. (2008) under Rosner’s model. They obtained 
Aq = 0.875 and Ai = 0.857. They also found that there is no evidence to reject the null 
hypothesis of equal cure rates. 
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Figure 3: Graph of the 81 empirical lengths of 90% HPD intervals generated from the 
case A = 0.0. 


Two striking features characterize this data, (i) The total sample sizes in the Cefaclor 
and Amoxicillin groups are extremely small, (ii) Table 2 is sparse with a zero cell 
in each group. As a result, the normal approximation used in classical analysis does 
not apply. In addition, the maximum likelihood estimates (MLEs) of Dallal’s reduced 
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g A=0.0 and m=25 

a 



A=0.5 and m=25 



I A=0.8 and ni=25 


f A=O.Oand m=10 

O 

E 



A=0.5 and m=10 



A=0.8 and m=10 




Figure 4: Graph of the 81 empirical lengths of 90% HPD intervals generated from the 
case A = 0.0,0.5,0.8. 


and saturated models sit on the boundary of the parameter space (leading to success 
probabilities not truly allowed by the product trinomial model). It is impossible to carry 
out frequentist inference or apply a likelihood-based model selection procedure. In this 
situation, a common ad-hoc adjustment is to add 1/2 to each cell count. A problem 
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Table 2: OME status after 14 days of Cefaclor or Amoxicillin treatments. 


Number of ears with 
OME being cured 

Treatment Group 

Amoxicillin 

Gefaclor 

0 

1 

0 

1 

0 

1 

2 

6 

3 

Total 

7 

4 


with this ad-hoc adjustment is that the total sample sizes are not integer numbers. The 
normal approximation still does not apply despite this adjustment. So the Bayesian 
methodology appears here to be one of the few alternatives. 


Table 3: MLEs under Dallal’s reduced and saturated models. 


Reduced Model 

No Adjustments 

U 

1 

V 

6/7 

7 

1/19 

0.95 

Ai 

57/70 

A 

-19/140 

R 

617 

iP 

3/13 


Id-hoc Adjustment 

u 

V 

7 

Ao 

Ai 

A 

R 

ijj 

10/11 

14/17 

1/11 

5/6 

77/102 

-4/51 

77/85 

lljVlh 


Saturated Model 

No Adjustments 

U 

V 

7o 

7i 

Ao 

Ai 

A 

R 

tp 

1 

6/7 

1/7 

0 

0.875 

6/7 

-1/56 

48/49 

6/7 

Ad-hoc Adjustment 

U 

10/11 

V 

14/17 

7o 

3/17 

7i 

1/27 

Ao 

17/22 

Ai 

27/34 

A 

4/187 

R 

297/289 

iP 

135/119 


Bayesian posterior estimates and 95% HPD intervals for Dallal’s model and the 
saturated model based on Jeffreys’, the reference, and the uniform priors and 100,000 
iterations are summarized in Tables 4 and 5. Under these three priors, the posterior 
means of U are far away from the parameter space boundaries as indicated by the non¬ 
zero median estimates (not provided for obvious reasons). The 95% HPD intervals for 
R under the uniform prior tend to be larger than those under Jeffreys’ prior or the 
reference prior. In addition, the uniform prior seems to perform the worst in terms of 
the deviance information criterion (DIG) while Jeffreys’ prior and the reference prior 
seem to perform similarly. The risk difference is essentially zero {BF^ = 1.965 and 
P(A > 0\D) = 0.550; BF^ = 1.818 and P{A > 0|Z1) = 0.541) and the two cure rates 
themselves are very high (above 75%). We retain the saturated model over Dallal’s 
reduced model according to the DIG and P( 7 i — 7 o < 0|iJ). However, the Bayes factors 
indicate minimal evidence against the null hypothesis Hq : 7 i = 7 o (PP^ = 0.682 and 
BF^ = 1.052). The 95% HPD interval for 71 — 70 also covers zero. According to the 
retained model, the conditional posterior probabilities of the cure rate at one site given 
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Table 4: Posterior estimates and credible intervals based on Dallal’s model. 


Jeffreys’ Prior 

P(A < 0\D) = 0.739, DIC=9.843, pD=1.613 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

U 

0.904 

0.119 

(0.649, 1.000) 

Al 

0.764 

0.125 

(0.515, 0.970) 

V 

0.819 

0.127 

(0.570, 1.000) 

A 

-0.079 

0.162 

(-0.417, 0.265) 

1 

0.076 

0.061 

(0.000, 0.198) 

R 

0.931 

0.245 

(0.468, 1.370) 

Ao 

0.842 

0.120 

(0.597, 0.998) 


1.013 

2.052 

(0.002, 3.423) 




Bernardo’s Prior 





P(A < OID) = 0.737, DIC= 

=9.953, pD=1.631 



mean 

std 

95% HPD 


mean 

std 

95% HPD 

u 

0.899 

0.124 

(0.635, 1.000) 

Al 

0.756 

0.128 

(0.505, 0.969) 

V 

0.812 

0.130 

(0.559, 0.999) 

A 

-0.082 

0.168 

(-0.439, 0.261) 

1 

0.077 

0.062 

(0.000, 0.199) 

R 

0.929 

0.260 

(0.469, 1.395) 

Ao 

0.838 

0.124 

(0.583, 0.997) 

tp 

1.030 

2.166 

(0.001, 3.486) 




Uniform Prior 





P(A < 0|£i) = 0.640, DIC= 

=10.378, pD=1.316 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

u 

0.833 

0.141 

(0.549, 1.000) 

Al 

0.713 

0.127 

(0.461, 0.935) 

V 

0.778 

0.131 

(0.524, 0.989) 

A 

-0.051 

0.177 

(-0.400, 0.317) 

1 

0.095 

0.066 

(0.003, 0.224) 

R 

0.973 

0.303 

(0.453, 1.536) 

Ao 

0.764 

0.136 

(0.491, 0.977) 

Ip 

1.200 

1.766 

(0.016, 3.796) 


the other site was cured are very high and slightly higher in the Amoxicillin group, 
therefore this dependency cannot be ignored. There is also a noticeable discrepancy 
between the correlation coefficients for the Zijk variables for the two treatment groups, 


5.2 Bayesian Analysis of Bilateral Data with Large Sample Data 

Postlethwaite et al. (2008) considered a two-arm multi-centre double-blind randomized 
trial where 168 diffuse scleroderma patients are randomized to one of two groups to 
receive either oral native collagen at a dose of 500g/day or a similar appearing placebo. 
The total duration of the treatment phase was 12 months with an additional visit at 
month 15 for safety follow-up. Rheumatologists routinely examine both the left and right 
feet, forearms, hands, fingers, legs, thighs, and upper arms of the patient and assign a 
modified Rodnan Skin Score (MRSS) score between 0 and 3, that is 0 for normal, 1 
for mild, 2 for moderate and 3 for severe skin thickening. The patient’s improvement 
at each body part level is recorded. After consultation with rheumatologists, a patient 
has improved at a body part level if the MRSS at the end of the trial is either zero or 
has dropped by two units or more from baseline. The goal is to test whether there is 
a significant difference in the improvement rates between the two groups at each body 
part level. Table 6 reports an examplary set of results from the trial for the forearms. 
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Table 5: Posterior estimates and credible intervals based on the saturated model. 


Jeffreys’ Prior 

P(A > 0|D) = 0.550, P( 7 i - 70 < 0|D) = 0.890, DIC=8.680, pD=1.6I9 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

U 

0.910 

0.112 

(0.670, 1.000) 

A 

0.021 

0.179 

(-0.342, 0.371) 

V 

0.824 

0.123 

(0.583, 0.998) 

R 

1.065 

0.299 

(0.548, 1.628) 

70 

0.192 

0.145 

(0.000, 0.481) 

V’ 

2.633 

7.094 

(0.001, 8.849) 

7i 

0.039 

0.055 

(0.000, 0.151) 

S 

-0.132 

0.167 

(-0.491, 0.200) 

Ao 

0.773 

0.128 

(0.532, 0.988) 

A 7 

-0.153 

0.155 

(-0.495, 0.114) 

Ai 

0.795 

0.125 

(0.553, 0.993) 








Bernardo’s Prior 



P(A > 0|D) = 0.541, P( 7 i - 70 < 0|D) 

= 0.889, DIC=8.858, pD=1.645 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

U 

0.900 

0.122 

(0.638, 1.000) 

A 

0.019 

0.186 

(-0.361, 0.393) 

V 

0.813 

0.129 

(0.560, 0.999) 

R 

1.067 

0.326 

(0.508, 1.667) 

70 

0.193 

0.145 

(0.000, 0.484) 

tp 

2.574 

6.512 

(0.004, 8.772) 

71 

0.040 

0.055 

(0.000, 0.153) 

S 

-0.135 

0.176 

(-0.502, 0.221) 

Ao 

0.765 

0.134 

(0.508, 0.988) 

Aj 

-0.153 

0.155 

(-0.506, 0.106) 

Ai 

0.784 

0.130 

(0.531, 0.989) 








Uniform Prior 



P(A > 0|D) = 0.583, P( 7 i - 70 < 0|D) = 

0.854, DIC=10.119, pD=1.352 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

U 

0.833 

0.141 

(0.548, 1.000) 

A 

0.040 

0.192 

(-0.339, 0.418) 

V 

0.777 

0.132 

(0.522, 0.989) 

R 

1.116 

0.386 

(0.476, 1.831) 

70 

0.231 

0.148 

(0.007, 0.521) 

V’ 

2.161 

3.727 

(0.022, 6.814) 

7i 

0.075 

0.073 

(0.000, 0.224) 


-0.110 

0.186 

(-0.484, 0.267) 

Ao 

0.686 

0.140 

(0.414, 0.938) 

A 7 

-0.157 

0.165 

(-0.515, 0.142) 

Ai 

0.725 

0.131 

(0.470, 0.956) 






In blocks 1 and 2 of Table 7, we provide Bayesian posterior estimates and cred¬ 
ible intervals for Dallal’s model based on Jeffreys’ and Bernardo’s priors. The re¬ 
sults are based on 100,000 posterior simulations. In block 3 of Table 7, we provide 
frequentist MLEs, Wald FCIs, Aikaike information criterion (AIC) and Bayesian in¬ 
formation criterion (BIC), along with large-sample f^st for Hq : Xq = Ai and 

Hq : 7 i = 7 o. For this data, the total sample size in each group is moderate. Al¬ 
though four of the cells have counts less than 5, the tallies that matter here are the 
subtotals, mi+ = 7 ,m 2 + = 9,mio + rn 2 o = 6 , moo = 55, mu -I-m 2 i = 10,moi = 36, 
which are all greater than 5. We provide only the results for Dallal’s reduced model. 
Indeed, Bayes factors indicate the null hypothesis, Hq : 71 = 70 , is supported. In this 
example, the treatment cure rate has decreased the disease risk by two fold. The risk 
of disease in both treatment and control groups remain high (above 80%). The null 
hypothesis of equality of the treatment cure rates, Hq : Ao = Ai is supported by our 
analysis. Overall Dallal’s model and the saturated model give similar results and there 
is no difference between using Jeffreys’ prior (e.g., R = 2.479 and BCI: (0.640, 5.007)) 
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Table 6: Number of scleroderma patients whose forearm MRSS decreased by 2 or 3, or 
has 0 MRSS at month 15. 


Number of forearms 
with improvement 

Treatment Group 

Collagen 

Placebo 

0 

36 

55 

1 

4 

3 

2 

6 

3 

Total 

46 

61 


Table 7: Posterior estimates and credible intervals as well as frequentist MLEs and 
conhdence intervals based on Dallal’s model. 


P{A > 0\D) 

= 0.956 

Jeffreys’ Prior 

, DIC=18.563, pD=2.878, BF^ = 1.607, BF^ = 3.733 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

U 

0.107 

0.039 

(0.037, 0.186) 

Al 

0.178 

0.049 

(0.087, 0.277) 

V 

0.229 

0.061 

(0.114, 0.348) 

A 

0.095 

0.057 

(-0.017, 0.208) 

7 

0.290 

0.099 

(0.109, 0.486) 

R 

2.479 

1.338 

(0.640, 5.007) 

^0 

0.084 

0.031 

(0.028, 0.147) 

tp 

2.855 

1.746 

(0.566, 6.124) 




Bernardo’s Prior 



P(A > 0|£i) 

= 0.957, DIC=18.523, pD=2.882, REa = 1.526, BF^ = 2.368 


mean 

std 

95% HPD 


mean 

std 

95% HPD 

P 

0.104 

0.038 

(0.037, 0.182) 

Al 

0.174 

0.049 

(0.085, 0.272) 

V 

0.223 

0.060 

(0.111, 0.341) 

A 

0.092 

0.056 

(-0.015, 0.204) 

7 

0.291 

0.099 

(0.110, 0.487) 

R 

2.481 

1.318 

(0.630, 5.004) 

^0 

0.081 

0.031 

(0.027, 0.143) 

tjj 

2.846 

1.707 

(0.582, 6.115) 




Likelihood 




AIC 

=18.712, BIC=26.730, pD=3 

^ xl = 

2.897, 

= 0.152 


mle 

std 

95% Cl 


mle 

std 

95% Cl 

u 

0.098 

0.038 

(0.024, 0.173) 

Al 

0.170 

0.049 

(0.073, 0.267) 

V 

0.217 

0.061 

(0.098, 0.337) 

A 

0.093 

0.057 

(-0.019, 0.205) 

7 

0.280 

0.102 

(0.081, 0.479) 

R 

2.210 

1.057 

(0.866, 5.641) 

Ao 

0.077 

0.030 

(0.017, 0.136) 

V' 

2.458 

1.323 

(0.855, 7.062) 


or Bernardo’s prior (e.g., R = 2.481 and BCI: (0.630, 5.004)), although results from 
Bernardo’s reference prior are easier to compute. Note that our posterior estimates and 
HPD intervals for the risk difference are in line with the results in Pei et al. (2010) 
obtained under the equal correlation model (A = 0.0970, Cl: (—0.0214,0.2217)). 

6 General Classes of Prior Distributions 

The reference priors and the uniform prior discussed earlier can be embedded in the 
family of prior distributions 
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2“ Ag»-^[l-(l + 7)Ao]'° ^ 


B(a, /3) (1 + ry)a+P-ao-ai 


Ar-^ [1-(1+7 )Ai]^-^ 
B(ai,foi) 


B(ao, &o) 


( 7 , Ao, Ai) G n, 


( 12 ) 


which is equivalent to stating that 


2 “ 7“-1(1-7)/5-1m“ 0-1 )hi-l 

= --57-TT-• (13) 


B(a,/3) (1+7)“+/3 B(ao,M 


B(ai, 61 ) 


Another representation of this class of prior distributions is through the following hi- 

2“ 7““^ (1-7)^“^ I-7 

erarchical model: (a) 7 - f{y) = — —— ————— or equivalently —- 

B(Q:,d) (l + 7 )“+P 1+7 

Be(/ 3 , a), (b) A0I7 -- Be ( ao, 6 o; 0 , —— ), and (c) A1I7 -- Be ( ai, 61; 0 , —— ), where 

V 1 + 7 / V 1 + 7 / 

Be(Q:, / 3 ; /, u) stands for a Beta random variable with shape parameters a and /3 defined 
on the interval {l,u). 

Jeffreys’ prior distribution for the parameterization (7, U, V) suggests a larger family 
of conjugate prior distributions for (7, [/, y), namely, 

B(a,/ 3 ) (1 + 7)“+^ ^ 

with 0 < 7, M, u < 1 and a, /?, oq, &01017 ^1 > 0, which translates into the prior distribu¬ 
tion 


7 r(7, Ao, Ai) oc 


+£+(Ao + .A.)V^Arll-(l+7)Aol‘”- 


B(a,/3) ( 1 -^7)“+/3-“o-oi-d 

xAr-'[l-(l+7)Ai]'^-' 


(7) Ao, Ai) G fl. 


(14) 


7 Concluding Remarks 

V -u 

Using the parameterization ( 7 , [/, U), it can be deduced that A = -. This result 

1 + 7 

highlights a direct dependence of the risk difference on the nuisance parameter 7 . This 
result also points out the main difference between the risk difference in 3 x 2 bilateral 
data and the risk difference, U — [/, in an ordinary 2 x 2 table where one collects a 
single measurement per subject. In other words, the divisor 1 + 7 is the term that 
connects the two disease risks in the bilateral data context. Indeed, according to the 
expression of the likelihood, +( 7 , t/, U), given in Appendix I, the parameters U and V 
can be interpreted as the proportions of cases with one or more body part(s) cured in 
the placebo and treatment groups. On the opposite side, Aq and Ai are interpreted as 
the proportions of body parts cured in the placebo and treatment groups. So unlike U 
and V, one individual can contribute twice in the computation of Aq and Ai. 
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Table 8: Relevant summary statistics for the risk ratio. 


Group 


Numbers of 


cured organs Treatment 


Control 


0 

1 or 2 


moi (1 - U) 
m+i - mpi {U) 


moo (1 - V) 
m+o - mpo (V) 


Another beneht of the parameterization (7, U, V) is that it shows that the risk 



frequentist and Bayesian inferences do not depend on 7 and are easier to compute. 
Therefore, the risk ratio has a technical advantage on the risk difference. Moreover, the 
definition of i? in 3 x 2 bilateral data coincides with the definition R from the 2 x 2 
binary table in Table 8. Actually, the 3 x 2 bilateral table can be replaced by Table 8 
when the focus is on R. For these reasons, our choice of the parameter in a bilateral 
data design is the risk ratio. 

While there are numerous frequentist papers dealing with bilateral data, this work 
remains incomplete. For example, Dallal’s model used in this paper has not been in¬ 
vestigated in the frequentist literature although other models have and there are no 
Bayesian treatments of the problem. Although the risk ratio and the odds ratio are well 
established parameters in medical settings, they do not appear in the bilateral data 
literature. A clear advantage of the risk ratio over the commonly used risk difference 
is that inference does not involve the nuisance parameter 7. In addition, the risk ratio 
remains unchanged when going from an ordinary 2x2 binary table where all observa¬ 
tions are independent to a 3 x 2 binary table where observations taken from the same 
subjects are correlated. The risk ratio and the odds ratio open the door for bilateral 
data to be studied under a prospective scheme or a retrospective scheme using an ob¬ 
servational study paradigm. The presentation exposed here has taken into account all 
these inconveniences and has provided a broad discussion of the Bayesian framework 
both from the point of view of the tests of hypotheses as well as the estimation of the 
key model parameters. We have added a simulation study to empirically compare the 
effectiveness of Bayesian methods against themselves as well as frequentist methods in 
the context of small, moderate and large sample sizes and for a wide range of A values. 
For example, we have found that Jeffreys’ prior and the reference prior tend to perform 
similarly. We have also found that frequentist methods tend to perform very poorly 
when A is small and the sample size is small or moderate. The uniform prior has the 
best overall property when it comes to estimating the parameter 7. We have concluded 
our work with two detailed case studies, one of which shows that it is impossible to carry 
out frequentist inference given that some of the parameters sit on the boundary of the 
parameter space or the normal approximation is not accurate. Our Bayesian framework 
works remarkably well in these situations as well as the large sample cases. 

When subject level bilateral data with covariates are available, Dallal’s regression 
model can be developed in order to incorporate covariates. As discussed in Section 3 . 2 , 
7, [/ = (1 -|- 7)Ao, and V = {1 + 7)Ai are unconstrained and the parameter space is 
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(0,1) for each of these three transformed parameters. Therefore, a logistic regression 
model can be assumed for each of 7 , U, and V. Consequently, Jeffreys’ prior and the 
reference prior can be derived. However, the computational and theoretical properties 
of these priors need to be carefully examined. The development of Dallal’s regression 
model deserves a future research project, which is currently under investigation. 
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Appendix 1: Derivation of Jeffreys’ Prior 

Set {7 = (1 + 7 )Ao and U = (1 + 7 )Ai. Under this new parametrization, the parameter 
space reduces to the interval [0,1] for each of the three parameters. Theorem A.l empha¬ 
sizes the role of a reparametrization technique in obtaining Jeffreys’ prior distribution 
in the original parametrization after having derived Jeffreys’ prior in the reparametrized 
space. 


Theorem A.l (Consistency Under Reparametrization). Consider a model 
dJl = {p{X\0), X € X,9 € 0} and let be an invertible transformation of 9. Then, the 
Jeffreys’ prior corresponding to the parameter (j), 7r{(j)), is that induced by the Jeffreys’ 
prior density of 9, Tr{9). We prove Proposition 3.1 and 3.2 below. 


Proof. As a function of ( 7 , U, V), the likelihood function simplifies to 


U, V) oc 


' ^ jjra-in+m^.n Q 

(X _(_ ^ ' '' ' 


This parametrization splits the likelihood function into three unrelated pieces, each 
piece related to a single parameter. This makes it very easy to derive Jeffreys’ prior. 
The first and second derivatives of the log-likelihood function with respect to 7 , U, and 
V are 


5 ^( 7 ) _ mi+ m2+ mi+ + m2+ d‘^l{y)_ mi+ TO2+ mi+J- 7712 + 
dy 7 1 — 7 1 + 7 ’ 7 ^ (1 — 7 )^ (1 + 7 )^ 

dl{U) _ {m2o + mio) moo dH{U) _ (m 2 o + mio) moo 

du ~ u i-u’ du^ ~ c 72 (1 - uy’ 

dl(y) _ (rn 2 i + mil) rnoi d‘^l{V)_ {m 2 i+mii) toq! 

dV ~ V 1 - V’ dV^ ~ U 2 (1 - U) 2 - 


Let r = 


m+i 

m+o 


be the ratio of the sample size in the two groups. Thus, we have 


r 52 /( 7 )! 

2 m+o (/7 + rV) 

52 /( 1 /) 1 

m+o 

57 ^ . 

7(1-7)(1 + 7)^’ 

5/72 

1 

I 


r d^KV)] 

. 5U2 . 


m+i 

U(l-U)- 
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Hence, Jeffreys’ prior under the parameterization ( 7 , U, V) is 
Trj{'y,u,v) oc 


(it + rv) 


7(1 — 7)(1 + 7 )^it(l — it)i;(l — v) ’ 


0 < 7 , It, !> < 1 , 


and it is proper. Theorem A.l implies that Jeffreys’ prior density in the parameterization 
( 7 , Ao, Ai) is 


7’‘j(7:-^o,Ai) oc 
where ( 7 , Ao,Ai) G ft. 


(1 +7)(Ao + rAi) 


7(1 - 7 )AoAi [1 - (1 + 7 )Ao] [1 - (1 + 7 )Ai] ’ 


□ 


Appendix 2: Derivation of Reference Priors 


Proof. We first derive the joint reference prior for the parameterization ( 7 , U, V) and 
then transform back to derive the joint reference prior for ( 7 , Aq, Ai). See Yang (1995) 
for justification. We also adopt the notation in Yang (1995). 


Case 1: Consider the group ordering {U,V} and then 7 or {V, t/} and then 7 . We have 

TO+o ^ _ m+i ^ _ 2m+Q{u + rv) 


hi = 


u(l — u) 


i;(l — v) 


and /13 = 


7(1-7)(1 + 7)2 


Thus, we obtain 


1 exp{i/log(/ii/i2)d7} 

7rii(MW,7) = ri, n/2, -Ti -^-Y’ 

^ ^ "*^exp|- JJJ \og{hih2)d'ydudv\ 

_ ^2 71 / 2 - 1(1 _ ^)l/ 2 -l ul/ 2 -l(l _ „l/ 2 -l(l _ ^)l/ 2 -l 

TT (I + 7 ) B(1/2,1/2) B(1/2,1/2) 


7r)j(7, Ao,Ai) 

= A^y/2-l(l_^)l/2- 

TT 


Ay"-'[l-(l+7)Ao]'^' 'a;/"-'[1-(1 + 7)Ai]'/" ' 


B(l/2,l/2) 


B(l/2,l/2) 


Case 2: Consider the group ordering of { 7 } and then {U, Y} or the ordering { 7 } and 

2 m+o(u + ri;) m+o m+i 

-jz -Tyr-— 77 , h 2 = — -r, hs = — --. Thus, the 

7(1 — 7)(1 + 7)2 it(l — u ) i ;(1 — u ) 

reverse reference prior is 


then {V,U}. We have hi = 


2 \hih2\^^'^ exp{lJJlog{h3)dudv} 

^ / |hih 2 |i/^du(ii; exp { i///log(h 3 )duc?i;d 7 } 

ul/2-l(l _ u)l/2-l ^l/2-l(l _ ^)l/2-l 72 71/2-1(1 - 7)1/2-1 
B(l/2,l/2) B(1/2,1/2) tt (I+ 7 ) 

=7r]j(u,i;,7). 
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Case 3: Here one starts with some subjective joint prior distribution for the parameters 
for which one has a good knowledge of and for the other parameters, one uses a non- 
informative prior distribution to reflect the lack of knowledge. Assume there is prior 
evidence for assuming the following conditional joint distribution 


^mao,a,i.)^ ■ o<Ao.A.< ' 


/ _ \ao+bo —1 / \ai+6i —1 • 


1 + 7 


about Ao and Ai given 7 . That is, 


A 0 I 7 ~ Be ^ao,&o;0, ^ ^ai,6i;0, 

and both Aq and Ai are conditionally independent given 7 . But one has no idea about a 
prior for 7 . A solution to this problem is to use a reference prior for 7 . Under our partial 
prior specification, the reference prior is Jeffreys’ prior associated to the integrated 
likelihood (integrating out U and V), 


L{y) oc 




(1 _|_ ^'jmi++m2+- ■ 

Hence, the reference prior corresponding to this integrated likelihood is 


3 ^ /w+oA(Ao) + to+iU(Ai) 




1 


OC 


7(1-7)(1 + 7)2 Y ^(1 _^)(i+^)2’ 

and it is proper. Thus, the joint prior distribution over H is then 

7r|(7, Ao, Ai) 

^y/2-i(i ^)i/2-i A“«-i - Ao) ” A^-^ (i^ - A] 


^ ( 1 ^)““^'“ 'B(ao,&o) ^ B(ai, 6 i) 

and it belongs to the family of prior distributions discussed in ( 12 ). 

Case 4: Similarly, in the second setup of partial prior specification, one assumes that 
there is available a family of prior distributions for 7 , for example. 


, ai+bi —1 


= 


2/A 

B(^,Ia) (l+yj/A+i- 


and one would like to find the joint reference prior distribution for the pair (Ao,Ai) 
conditional on 7 . The proposed reference prior under this partial prior specification is 


4 (Ao,Ai| 7 ) oc |det(5)r/" = ^+ 22^33 cx 


(1+7) 


AqAi [1 — (1 + 7 )Ao] [1 — (1 + 7 )Ai] 
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where S is the lower 2x2 left corner matrix of the Fisher information matrix. Hence, 
the joint prior distribution of ( 7 , Aq, Ai) over fl is 

7rfl(7,Ao,Ai) 

^ [1 - (1 + 7 ) Ao ]^^^~^ [1 - (1 + 

(1 + 7 )m+-i B(1/2,1/2) B(1/2,1/2) 

and it belongs again to the family of prior distributions discussed in ( 12 ). □ 


Appendix 3: Results for Posterior Calculation 

Proposition A.l. Let 7 have density fij) = — - r— - -r-^r -, 0 < 7 < 1. 


B(/r,z/) ( 1 + 7 )^'+^' 


pn4‘ 


Then, the density of (j) = logit{'j) + log(2) is f{(j)) = — -- 

known to be the density of logit{p), where p ~ Be(^, i/). 


, which is well 


O0-Iog(2) 


Proof 7 = 
dj 


. We have 1 — 7 = 


1 


1 _|_ g<ii-log(2) 


^ 1 + 2 e'^-'°s( 2 ) 

1 + 7 = j 

' 1 _(_ g0-log(2) 


. Thus, we obtain 


d(j) (1 + 

2^ g(Ai-l)( 0 -log( 2 )) 


/(</') = 


(X _|_ g0-log(2))M+i^ g0-log(2) 


B(xt, z/) (1 + (1 + (1 + (1 + 


21^ 


„/i(<^-log(2)) 




B(/r, v) (1 + 2 e‘^-i°g( 2 ))M+!^ (x + e‘^)''+‘' 


.□ 


2 M ^M-1 (X - 

Proposition A.2. Let 7 have density fipf) = — --——- 7 — -, 0 < 7 < 1. 


1 — "T 

Then, tt =-~ Be(z^, xt). 

1+7 


B(xt, z/) (1 + 7 ) 7 *+'" 


1 — TT 27r 2 d'y —2 

Proof. We have 7 = ——. So 1 -7 = ——, 1 +7 = ——, and — = — —Thus, 

l+TT l + TT l+TT dlT ( 1 + 

we have 


/(t^) = 


2 + (1 - tt)'"-! (1 + 7r)++^ 2 _ 1 

B(/i, z/) (1 + 71 ) 7 *“! (1 + 77 )“"“^ 27*+*" (1 + tt)^ B(z^,/ i) 


7r‘^-^(l-7r)'^-i 


□ 


The marginal prior distribution 71 ( 7 ) has the following properties. When p = y = 1, 
77(7) is decreasing. When X7 = 1 and v < 1, 77(7) is U-shaped, the anti-mode being at 
y. When /r = 1 and y > 1 or y = 1 and X7 < 1 or x* > 1 and y < 1, 77 ( 7 ) is decreasing. 
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When 1 / = 1 and 1 < // < 3, is unimodal and the mode is at 7 = ^ ^ ^ . When 

1 ^ = 1 and > 3, it is J-shaped. When ijl,v > 1, 7 r( 7 ) is unimodal and the mode is at 

7 = -where A = (2i/ + /r — 1)^ — 8(/x — 1). When ^ < 1 and < 1, 

2i/ + /r - 1 + V A 


7 r( 7 ) is U-shaped, and the anti-mode is at 7 


2 z/ -t- /r — 1 -j- "x/A 
4 


Appendix 4: Posterior Distribution of A and R 


To derive the posterior distribution of A and R, we consider the parameterization 
( 7 , U, V). Under this parameterization, the posterior distributions of interest belong to 
the family 

1 omi+S-o; /'i _ \m2++i3-l 

Tr(y,u,v\D) = 77:57-77-77-r-7- , ,0 — {u + rvY 

K B{mi+ + a,m2+ + P) (1-|-7)™i+-i-"*2+-i-a-i-/3 ^ ^ 

ymw+'m20+ao — l _ y^moo-|-&o —1 ymii+m2i+ai-l _ y'^moi+bi—l 

B(toio + m 20 + ao,moo-f 6 o) B(toii-|- m 2 i-I-oi, moi-b 6 i) 


where 0 < 7 ,m,z; < 1 and K is the normalizing constant. The choice d = 1/2, oo = 
Bq = ai = bi = 1/2 corresponds to Jeffreys’ posterior distribution and the choice 
d = 0, oo = 6 o = ai = bi = 1/2 corresponds to Bernardo’s posterior distribution. 

Proposition A. 3. The posterior distribution of the risk difference, A = Ai — Aq = 
V -U 

-, has the complex integral form 

1-1-7 


r(A|iA) 


1 

K 


1 -|A| 

|A| 


2 mi + -|-a 


r^m-i++a-l _ ,^'^m2++h- 


I-HAI B(toi+ -b a, 7712+ -b /3) (1 -b 7 )™i++"*2++“+/3 


—(1+7) A) 
/ max( 0 , —(1+7) A) 


-b (1 -b 7)A^ 


/ ~.dipnio+m2o+ao-l(l_y)'moo+bo-i 

1^(1 -b r)u + r(l -b 7)Aj 

( 


mii+m 2 i+ai-l 


1 — U — 


B(mio + TO 20 + oo, moo + bo) 

moi+bi-l 


(1 +7)A) 


B(mii -b m2i -b oi, moi -b 61) 


■ du dy. 


lUe also have 


P[Xi — Xq > Ao|I?) 


2 "*i + + 5 7 7 7 7 "*i+ + 5-1(1 - y)™ 2 + + 5-l 

B(mi+ -b i,m 2 + -b i) Jo Jo Jv>u+Ao/{l+'y) (1 -b7)'"i + +™ 2 + + l 

_|_ rv)‘^ yyi-lQ+m20 + ^-l _ y^moo-l-|-l ymii+m2i + ^-l _ ^'^moi + ^-1 

K B(mio-b m 2 o-b 5 ,moo + 5 ) B(mii-b m 2 i-b 5 ,moi-b 


dvdudy. 


When Aq = 0 , this expression depends no longer on 7 and when d = 0 it is even 
simpler. 
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Ai y 

Proposition A.4. The posterior distribution of the risk ratio, R = — = does not 
depend on 7 and it has the simpler integral form 


-k{R\D) = < 


j^mii+m2i+ai-l -(_ rR)‘^ 
KK* 

J^-{mio+m2o+ai+d+l) _j_ 

KK* 


u 


mi++m 2 ++d+ao+ai —1 


(1-w) 


moo+bo — i 


X (1 — Ru) 


moi+bi-l 


du, 


0 < i? < 1, 


^l++m2++d+ao+al-l _ ^^moi+&i-l 


X (1 - Rv) 


"1.00+&0—1 


dv , 


i? > 1, 


where K* = B(mio + m 2 o + oq, toqo + 6 o)B(mii + m 2 i + ai, moi + 61 ). When d = 0, it 
can be shown that is unimodal when moo bind toqi bire positive. We also have 


P{R > Ro\D) = 


'fa I rv ^2 ^(1 2 ^ 


0 Jv>Rqi 


K 


B(mio + TO 20 + 5 , "loo + \) 


mii+m 2 i + -l-l 


( 1 - 1 ^) 


moi + i —1 


B(mii + 77121 + "loi + 5) 


dvdu. 










